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A model of mobile, charged ion channels in a fluid membrane is studied. The channels 
may switch between an open and a closed state according to a simple two-state kinetics 
with constant rates. The effective electrophoretic charge and the diffusion constant 
of the channels may be different in the closed and in the open state. The system is 
modeled by densities of channel species, obeying simple equations of electro-diffusion. 
The lateral transmembrane voltage profile is determined from a cable-type equation. 
Bifurcations from the homogeneous, stationary state appear as hard-mode, soft-mode 
or hard-mode oscillatory transitions within physiologically reasonable ranges of model 
parameters. We study the dynamics beyond linear stability analysis and derive non- 
linear evolution equations near the transitions to stationary patterns. 



I. INTRODUCTION 

Spontaneous pattern formation of ion channels in cell 
membranes has been a subject of continuous interest dur- 
ing recent years Jl|-§| ■ Spatially modulated distributions 
of ion channels or pumps are ubiquitous in many cells. 
They are closely related to important biological functions 
like some early stages of morphogenesis (J] or the stor- 
age and processing of information in neural tissue ||. 
Many physiological regulation processes involve changes 
of ion fluxes through cell membranes on time scales rang- 
ing from milliseconds to hours ||. 

There have been attempts to explain the aggregation 
of channel proteins in a fluid membrane on the basis 
of thermodynamic models JlJ , [Q and references therein. 
Typically the system of channels is driven out of equilib- 
rium due to ionic concentration gradients and transmem- 
brane fluxes. Therefore, a number of authors have put 
forward models for spontaneous pattern formation based 
upon systems of semi-phenomenological, non-equilibrium 
equations of motion for the densities of channel proteins, 
the charge densities of ions and the intra- and extracellu- 
lar voltage Jj],[|J§ • In many important systems like axons 
or cells in neural or muscular tissue, there is a bulk vol- 
ume of aqueous solution of ions on one side of the mem- 
brane whereas on the other side, there is only a thin layer 
of electrolyte, whose thickness usually is in the range of 
10 • • • 100 nm. See Fig|| For a schematic picture of the sit- 
uation. Under these conditions, the models may be sim- 
plified further. Ion densities may be eliminated approxi- 
mately, leading to a coupling of the densities of channels 
to the lateral transmembrane voltage profile within that 
layer. The voltage profile has to be determined from a 
quasi one-dimensional or two-dimensional cable equation 
P,p[JlO|] . This contains transmembrane currents and thus 
channels are effectively coupled via the transmembrane 
voltage. 



In the present work we will follow this approach and 
analyze a model of charged ion channels, which may 
switch between a closed and an open state, thereby also 
changing their effective electrophoretic charge and their 
diffusion constant. The effective electrophoretic charge 
can be of both signs - although the proteins are usually 
only negatively charged - if the electro-osmotic effect is 
taken into account S. 

We consider this model as a natural, minimal model of 
mobile channels with kinetics of internal states. As the 
conformations, which change the states of an ion channel 
prLfni, may involve charge transport across the mem- 
brane Jl^,[l3| as well as binding to other molecules, for 
instance phosphorylation via protein kinase or immobi- 
lization due to binding to the cytoskeleton jT^|l4|], the 
state dependence of charge and diffusion constant should 
be taken into account. 

In (|] a special case of the present model has been stud- 
ied already. The authors of this work considered a reser- 
voir of closed channels in the well-stirred approximation 
without coupling to the lateral voltage profile. We will 
see below that our model leads to qualitatively different 
behavior in a wide range of parameters but still contains 
the results of ^ . The necessary limit is that closed chan- 
nels diffuse faster than the open ones and that their elec- 
trophoretic charge vanishes. The present model is also 
related to [|l5| , where a mixture of two different channels 
has been considered without a kinetic, which allows for 
transitions between the different species. 

The basic mechanism, which drives spontaneous pat- 
tern formation is very simple in models, which only in- 
volve mobile, charged ion channels and the lateral profile 
of the transmembrane voltage || . It is based on a feed- 
back loop during which the channel proteins drift due to 
lateral voltage gradients and at the same time modify 
the voltage gradients by Ohmic voltage drops, caused by 
transmembrane currents. Therefore, for example, neg- 
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atively charged channels, which cause currents directed 
into the thin layer of electrolyte have a tendency to aggre- 
gate. Diffusion will counteract this tendency and thence 
a variety of bifurcations of the homogeneous state to an 
inhomogeneous state becomes possible. We will show be- 
low that our model implies a number of different scenarios 
for the first bifurcation from the homogeneous state. 

Depending on parameter ranges either soft-mode in- 
stabilities arise, where periodic patterns appear with a 
wavenumber k c , which approaches zero at the transition. 
Or there will be different kinds of hard-mode instabili- 
ties, where patterns with k c ^ at the transition and 
spatio-temporal patterns emerge with non-vanishing k c 
and non- vanishing frequency f2 c . 

After introducing the model in the next section, we 
will give results of the linear stability analysis together 
with results of non-linear behavior close to transitions 
into stationary patterns in subsequent sections. For soft- 
mode instabilities, we derive a Cahn-Hilliard type equa- 
tion of motion for the slow modes near the transition. 
For hard-mode instabilities, an amplitude equation is de- 
rived, which allows us to separate parameter regions of 
forward and backward bifurcation or continuous and dis- 
continuous transitions respectively. 



II. THE MODEL 

We consider ion channels moving within a fluid mem- 
brane of size L x x L y which separates a thin layer of 
electrolyte of width d from an electrolytic bath, as it is 
sketched in Fig.[l| Positions r within the membrane are 
described by internal rectilinear coordinates r = (x,y), 
< x < L x , < y < L y . If we choose periodic bound- 
ary conditions in the y coordinate, the membrane corre- 
sponds to a cylindrical cable and is reminiscent of axonic 
or dendritic structures for L x > > L y . Note that we treat 
L y and d as independent quantities, so that the actual 
structure in cylindrical geometry looks like on the right 
hand side of Fig.0, with a sub-membrane layer of width d 
and a 'core' of radius (L y /2ir) — d, which is decoupled 
from the sub-membrane layer. For d = L y /2w the sub- 
membrane layer fills the entire interior of the cable. In 
the following we will not distinguish between these dif- 
ferent geometries, as the present work only deals with 
infinitely extended systems and hence effects of specific 
boundary conditions need not be taken into account. 

Ion channels may switch between an open state (V) 
and a closed state ('c') according to a simple mono- 
molecular chemical reaction scheme 
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with rates 7_ and 7+ . We are interested in a wide range 
of kinetic rates, from 10 3 /sec for voltage or ligand gated 
channels down to less than one per hour as it is observed 



for hormone regulated phosphorylation and dephospho- 
rylation of channel proteins. The distribution of chan- 
nels is described by smooth densities n r (r, t), r e {o, c}, 
which obey the equations 

d t n r (r, t) + Vj r (r, t) = a r [y+ri c (r, t) - 7_n (r, t)] (2) 

with o r — +1 for r =0 and oy = —1 for r =c. The rhs 
of this equation takes into account the reaction kinetics 
of ([!]). Current densities j r of the channels are assumed 
to be of the Nernst-Planck form 

j r (r,f)= -D r {Vn r (r,t) - 0n r (r,t) -g r E(r,t)} (3) 

where D r and q r refer to the constants of lateral diffusion 
and to the effective electrophoretic charges of open and 
closed channels, respectively. The coefficient is the 
inverse temperature and E(r) = — VV(r) denotes the 
lateral electric field. Note that we explicitly allow for 
different diffusion constants and different electrophoretic 
charges in the states and c. The simple form of the 
current densities given in (^) may be justified for rigid 
protein structures with charged protuberances extending 
in the intra- and extracellular electrolyte. As discussed 
in H) , the effective electrophoretic charges may be of both 
signs due to electro-osmotic effects, although proteins are 
usually negatively charged. 



As has been discussed in the literature 10 1, the Kelvin 
cable equation or its two dimensional analogue may be 
used to calculate the lateral variations of the transmem- 
brane potential provided that spatial variations of ion 
concentrations have negligible effects and that the char- 
acteristic length scales of lateral patterns are large com- 
pared to the width of the thin layer of electrolyte. In [jlO| , 
the cable equation is derived by integrating the Nernst- 
Planck equation of ion densities together with the Pois- 
son equation using appropriate boundary conditions at 
the membrane. In the following we use the 2-dimensional 
cable equation S 



C m d t V(r,t) = —V 2 V(r,t)-GV(r,t) 

Pe 

-\n (r,t){V(r,t)-E}. 



(4) 



C m denotes membrane capacitance per area, p e the resis- 
tivity of the electrolyte within the thin layer of width d, 
and Xn Q is the conductance of open channels. We have in- 
cluded a passive, homogeneously distributed transmem- 
brane conductance G. E is the reversal potential of ion 
fluxes through the open channels, which drives the sys- 
tem out of equilibrium. 

We are interested in the stability of the homogeneous 
and stationary solution V,n ,fi c of (0) 
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We consider the total density of channels n as con- 
stant and define T := 7+ + 7_. Furthermore, we in- 
troduce appropriately normalized deviations from the 
stationary and homogeneous solution as the new fields 
p = {n — n )/n , ( = (n c — n c )/n c for the open and 
closed channels, respectively, and $ = (iq {V — V)/V for 
the membrane potential. Next, we change to dimension- 
less lengths and times by using as units the typical decay 
length t v = [(p e /d)(Xfi + G)]~ 1/2 of V and a typical 
diffusion time r = i v / D a over this length. Then (§), (g) 
and (^) become dimensionless 



dtQ = V {Vg + (g + 1) V$} - 7 _(f> - C) 



(6) 



a t C = ^V{VC + g(C + l)V$} + 7+ (e-C), (?) 



e <9 t $ = [V 2 - l] $ - a<S>g - r)Q . 



(8) 



In these equations we use the ratio of diffusion constants 
D = D c /D , of electrophoretic charges q = q c /q , the 
relative conductivity of the membrane a — \h / (Xn +G) 
and a rescaled reversal potential rj = — a(l — a)f3q a E. 
The parameter e — C' m p e D /d compares the electrical 
relaxation (RC) time per area to the diffusion constant. 
Note that 7+ and 7_ have also been made dimensionless 
7± — > 7±r by multiplication with r. 

For cell membranes, the parameters in (g), (J3J) and (^) 
imply time-scales, which differ by several orders of mag- 
nitude. Typical parameter values are Jl2] C' m = 1/iF/m , 
p e ~ lflm, d ~ lOnm, and D r ~ 0.1(/im) 2 /sec for dif- 
fusion of mobile proteins. Channel conductances A are 
in the 10 _12 f2 _1 range, average channel densities are 
around 10 • • • 50 (fim) -2 and electrotonic lengths ly are 
of the order of a few /im. Then the unit r ~ 10 2 sec, 
whereas the electrical relaxation time is below 1 //sec and 
it is always the smallest time scale of the system, even 
compared to fast kinetics of voltage or ligand gated chan- 
nels with time constants in the msec range. The con- 
stant e ~ 10~ 7 . This separation of time scales justifies 
the use of the quasi-stationaryapproximation in (|8|) and 
hence we will put the rhs of (|Sj> equal to zero in the fol- 
lowing. For voltage or ligand gated channels, the time- 
scale set by the rate T is also very fast compared to r, 
r ~ 10 5 , whereas slower regulation processes have time- 
scales comparable to or larger than r. The limit T ^> 1 
will be referred to as fast reaction limit in the following. 

For reversal potentials \E\ ~ 0...100mV and elec- 
trophoretic charges of a few elementary units the param- 
eter rj is in the range \rj\ ~ • • • 10 under physiological 
conditions. 



III. LINEAR STABILITY ANALYSIS 



We linearize Q6|-|Sj) around the homogeneous and 
stationary solution (js)) and apply the plane wave 
Ansatz (g, £, $) = (gk, (k, $fc) exp{ik • x + uit} + c.c. 



As it was discussed above, within the quasi-stationary 
approximation potential fluctuations are proportional 
to fluctuations of open channels 



3>fe = - W(l + £: 2 ). 



(9) 



Inserting (^) into the linearized (^|) and (^) we get by 
abbreviating rjk — r//(l + k 2 ) the eigenvalue problem 



Qk 

Ck 



-7- 7- 
7+ -74 



k 2 l 1 - Vk 



gk 

-Dqi] k -D ) \ \ C, k 
The eigenvalues 

w± = -[P(k 2 ) ± v/Q(^)]/2 

p(k 2 ) = k 2 (i + D- Vk ) + r 

Q(k 2 ) = [k 2 (l-D- m ) +1 _- 1+ ] 
+ ADqj_T] k + 47_7 + 



(10) 



(11) 



determine linear stability. Let us note in passing that an 
extended model with voltage dependent reaction rates 
will lead to the same set of linearized equations. 

The solution (|5|) becomes unstable against small per- 
turbations if Rew + or Rew_ is positive. Note that os- 
cillatory unstable modes, which require Imw± 7^ and 
thus Q < 0, are possible in the parameter range qr] < 0. 

The decisive influence of the parameters rj and q on 
the qualitative behavior of solutions is already apparent 
in the fully non- linear system (|[g). The term rjg in (g) 
drives the curvature of the potential <i>, whereas q may 
change the relative direction of drift currents of the closed 
channels. For example, if open channels carry negative 
charge, their drift is directed towards depolarized regions 
of the membrane, where they accumulate. If a current 
through open channels is directed into the thin layer of 
electrolyte, this current will further depolarize the mem- 
brane and may cause an instability. If closed channels 
have a positive effective electrophoretic charge, the accu- 
mulated open channels will disperse again after closing. 
This may lead to an oscillatory behavior. By varying the 
relative signs of q and rj, different scenarios of aggrega- 
tion and dispersion of channels by lateral structures in $ 
may arise, which are coupled back into the (^) and may 
lead to stationary or oscillatory patterns. 

In the following we will hence discuss the four param- 
eter regions lZ ++ ={q > 0,r/ > 0}, 1Z ={q < 0, 77 < 0}, 

Tl + -={q > 0, 77 < 0} and Tl- + ={q < 0, rj > 0} sepa- 
rately. The quantity rj is treated as the primary control 
parameter. 
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A. Instabilities for qr\ > 

In these regions of parameters, Q(k 2 ) is always posi- 
tive, the rates 0J±(k 2 ) are real and thus plane wave so- 
lutions will either exponentially grow or decay in time. 
Instabilities signal the onset of stationary patterns. As 
cj+ > u)-, it is the rate w+, which will drive the solu- 
tion (||) unstable, if control parameters reach the stability 
boundary. Note that in the long wavelength limit k — > 0, 
uj + vanishes for all values of i],q,D,-f± like w+ = 0(k 2 ). 
The corresponding diffusion mode reflects conservation 
of the total number of channels. 

From ui+ui- = we get the neutrality condition 



m {k) = (i + k*) 



2 , D{k 2 + 1 -)+ 1+ 



D(fc 2 + <n-) + 7+ 



(12) 



The minimum of r]o(k) corresponds to the first onset of 
instabilities and the corresponding wavenumber k c is the 
first unstable mode (critical mode). A minimum of the 
neutral curve appears at k c ^ with 



k 2 = - g7 _ - 2± + v/ 7 _( g -l)(g 7 _+ 7+ /£>-l) (13) 



if the rhs of (|13j) is a positive, real number. By varying 
q, D or the rates 7 ±, the instability may switch between 
soft-mode (k c = 0) and hard-mode (k c ^ 0) type. 

Several statements are noticeable: (a) For D = 0, 
(immobile closed channels), as well as for q = 1, the 
neutrality condition, r/o = 1 + k 2 , is independent of 
the reaction rates and implies a soft mode instability 
at rj c = 1. A hard mode instability requires both mo- 
bile channels and different charges of closed and open 
states, (b) Hard mode instabilities are only possible for 
q < 1, i.e. q c < q a . (c) For D — > oo and q — we 
recover the result of (d) In the fast reaction limit 
( 7 _ >> 1,7+ >> 1,-D 7^ 0), only soft mode instabili- 
ties at rj c = (Z) 7 _/ 7+ + l)/(-D<7 7 _/ 7+ + 1) remain. Note 
that the critical control parameter stays small (0(1)) in 
the fast reaction limit, indicating that the instability may 
indeed be reached under physiological conditions. 

For reaction rates of order 1, typical curves of r)o(k) 
are shown in Fig.||. Whereas in 1Z++, both soft and 
hard mode instabilities may occur, a neutral curve in 

1Z does only appear for D\q\ > 7+/ 7 -, in which 

case r)o(k) is maximal at k = and exhibits poles at 
k 2 = —q-f- —j + /D. Thus, for electrophoretic charges q 
and q c of different signs, merely soft-mode instabilities 
at rj c < are possible as it is demonstrated in Fig.^.b). 
Note that only a finite band of unstable modes, limited 
to the interval (—k p ,k p ), will appear for all rj < 0. This 
unstable band grows like k p ~ T 1 / 2 in the fast reaction 
limit. Fig.|| shows the k 2 = lines in the £>-r)-plane for 
different values of q. For parameters (D,T) below each 
curve k 2 0. 

Setting Qk c = 1 and solving the linearized system for 
the eigenvector (gk e , (k c , &k a ) of the critical mode, we 
obtain Cfe c = (qk 2 + g 7 _ +j+/D)/{k 2 + 97 _ + 7+ /£>) and 



= -(k 2 c + 97- + l+/D)/[kl + qj- + 7+/£>). Taking 
a closer look at the signs we find sgn g^ a = sgn (fc c and 
sgngfc c = — sgn<E>fc c if q > 0. Hence for q > 0, spatial 
variations of open and closed channels are in phase but 
have a phase lag of tt with respect to the potential in 
the critical mode. But if q < there is no phase lag, as 
then sgn uik c — sgn Cfe c = sgn $fe c . In this case the critical 
wavenumber vanishes leading to Co — 1 and $o = ~ 7 h 
so that all spatial periods are in phase as expected from 
the qualitative discussion on the mechanism of pattern 
formation given in (O) after (0). 



B. Instabilities for qr) < 

The simplest behavior of the system appears in the 

regime lZ-\ For Imw± ^ 0, it is obvious from jll| ) 

that Reu± = — (l/2)P(fc 2 ) < 0. It is easy to show that 
linear stability also holds for Ima;± = and thus the 
homogeneous, stationary solution (^|) is linearly stable in 
this regime. 

It is only in 1Z |_, that oscillatory instabilities may 

arise, because Q{k 2 ) and P(k 2 ) can simultaneously take 
on negative values in this regime. If Q < 0, uj + = u*_ , and 
the growth rate of the perturbation is Rew± = —P(k 2 )/2. 
From the condition P(k 2 , rj) — we obtain the neutrality 
condition r}p(k) for oscillatory instabilities 



ri P (k) = (1 + D)(l + k z ) 1 + 



r 



(1 + D)k 2 

It attains its minimum at the critical wavenumber 



k P = 4 



r 



1 + D 

The critical value of the control parameter is 



(14) 



(15) 



(16) 



Note that P(k 2 ) does not depend on the ratio q of the 
electrophoretic charges and thus the neutral curve only 
depends on reaction rates and diffusion coefficients. The 
ratio q will, however, show up in the oscillation frequency 
f2 c = (l/2)-\/|<2(fc 2 ) at the critical point. 

The region 1Z |_ also contains parts with Q(k 2 ) > 0, 

where stationary unstable patterns may arise if the neu- 
tral curve determined from (^) is crossed. Hence the 
true critical point is located at 



Vc 



mm{rjo(k),r]p(k)} . 

k 



(17) 



The oscillatory instabilities will not be attainable in the 
fast reaction limit under physiological conditions because 
r)p c ~ r for large T. 

Fig.^ shows neutral curves in 1Z y together with the 

region Q(k 2 ) < for two cases, which differ by the value 
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of D. In Figjija, the unstable pattern is oscillating, in 
Fig.^j.b, it is stationary. FigJ| shows the dependence of 
the critical frequency tt c on D and T for different q. Note, 
that by decreasing D, the oscillations will set in with a 
finite frequency at the transition, whereas by increasing 
D the critical frequency f2 c starts out from f2 c = at the 
transition. 



IV. EFFECTS IN THE WEAKLY NONLINEAR 
REGIME 



We now turn to a discussion of the behavior of the sys- 
tem beyond linear stability analysis. This discussion is 
restricted to stationary patterns. A more detailed anal- 
ysis of oscillatory patterns will be given elsewhere. 

In the vicinity of the critical control parameter r\ c and 
the critical wavenumber k c the linear dispersion of the 
slow mode k 2 ) may be expanded to obtain the lin- 

ear part of a non-linear evolution equation for the slow 
mode near the bifurcation point. From (O) one gets 



+ £ (k - k c f 



-0[(k-k c ) 4 



(18) 



where z a = 4 for soft-mode instabilities (a = s) and z a — 
2 for hard-mode instabilities (a = h). The expressions 
for time scales r a and length scales £ a can be found in 
Appendix A for soft-mode and in Appendix B for hard- 
mode instabilities. 



A. Soft-Mode Instabilities 

To proceed beyond the linearized dynamics near a soft- 
mode instability, it is convenient to introduce the variable 
u = (7 + g + 7_e)/r, which corresponds to fluctuations of 
the total channel density, u — (n Q + n c — n) /n, and which 
contains the critical slow mode as rj — > r) c . Furthermore 
we introduce v — (g — C)/I\ which remains fast relative 
to u near the bifurcation. We get an equation for the 
slowly varying part of u(x, t) from a gradient expansion, 
which is outlined in Appendix A. 

If we take into account all terms of 0(V 2 ) and include 
the 0(V 4 ) terms, which arise from the linear dispersion 
(|l8|), the resulting equation for u takes on the form of a 
dynamical Cahn-Hilliard ( or 'model B') equation 



d t u 



-i v s 



SF{u) 



+ 0(V V) , 



(19) 



with an effective potential F(u) — f[f(u) + 

(e s /2)(\7u) 2 ]d 2 r. 

The local part f(u) is given by 



ft \ ^ cU , ^ 



(2 - a)u 



(20) 



1 u) ln(l + cm) 

a 



Equation ([19|) has to be solved under the constraint of 
conservation of channel number, J u d 2 r — 0, which may 
be taken into account by an appropriate Lagrange mul- 
tiplier term — Au in /(it). Note that the dynamics only 
depends on the ratio r] c /rj and on the parameter a and 
thus exhibits a particularly simple universality for all per- 
missible values of D, g, r y±. 

Fig|i| shows the effective potentials for a = 0.6 close to 
the transition. rj c becomes the spinodal. The discussion 
of the dynamical behavior, including typical coarsening 
in the case of quenches can be found in standard text- 
books, see for instance fl6[| . 



B. Hard-Mode Instabilities 

In order to obtain an analytic description of the weakly 
nonlinear regime for hard-mode instabilities, we use a 
standard multi-scale perturbative approach, which re- 
sults in an amplitude equation JlTj . The approach starts 
from the Ansatz 



u, 



[ X (X,Y,T)e^ x + X *(X,Y 1 T)e- 



(21) 



for the fluctuations u = (g, £, $). The amplitude \ de- 
pends upon the scaled variables X = ax, Y 



and T = a 2 t. The scale a becomes small, if the control 
parameter 77 approaches its critical value r\ c from above, 
a 2 = (Tj—ricJ/rjc. For convenience we sketch the procedure 
in Appendix B. For one-dimensional, stationary patterns 
the perturbation theory gives for the amplitude \ in the 
supercritical vicinity of i] c the equation 



ThdtX = 



Vc 



X-J\X?X- (22) 



The lengthy expression, which gives J as a function of 
the model parameters, is given in the Appendix B. Fig.|v| 
displays the sign of J in the D-F-plane for different val- 
ues of a. Figjy] shows the distribution of the sign of J 
within the _D-r-plane for different values of q and a. Pos- 
itive J corresponds to a continuous transition (forward 
bifurcation) negative J to a discontinuous one (backward 
bifurcation) . 

Note that for a < 0.5, Fig.|y|.a and Fig |y|.b, there is 
inside of the domain of k 2 > only one region of pos- 
itive J, i.e. continuous transitions, which is gcncrically 
bounded by the domain of negative J, i.e. discontinuous 
transitions. In the other regime a > 0.5, Fig.^.c and 
Fig ^.d, the system shows a reentrant behavior, as then 
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an additional domain with J > appears below the do- 
main of J < 0. This means, coming from large values of 
D one leaves the area of continuous transitions, crosses 
the domain of discontinuous transitions and again enters 
a domain with J > 0, before finally the k 2 — 0-line is 
crossed. Parallel to the T axis the generic behavior is 
similar. 



ter regions of forward and of backward bifurcations may 
be distinguished by the sign of the non-linear coupling 
parameter. 

The presented model includes and extends several pre- 
viously studied models and may be considered as a min- 
imal model for charged ion channels with internal state 
kinetics in a fluid membrane. 



V. CONCLUSIONS 

We have studied a simple model of spontaneous pat- 
tern formation of mobile ions in a fluid membrane. The 
channel proteins may switch between an open and a 
closed state according to a simple two state reaction ki- 
netics with constant rates. The effective electrophoretic 
charge and the diffusion constant of the channel proteins 
are assumed to be state dependent. The reversal poten- 
tial E of ions, which may pass the open channels consti- 
tutes a non-equilibrium driving force. The characteriza- 
tion of a particular model is completed by the additional 
parameters D = D c /D a and q = q c /q , the reaction rates 
7± and the fraction a of conductance due to the aver- 
age number of open channels to the average total trans- 
membrane conductivity. Varying these parameters leads 
to a number of distinct scenarios for first bifurcations 
from the stationary and homogeneous state. Depending 
on the sign of the control parameter 77 = a(a — l)(3q E 
and the ratio q — q c /q of electrophoretic charges, four 
qualitatively different regions of parameters have been 
identified from linear stability analysis: (a) for q > 
and 77 > (1Z++) soft- or hard mode instabilities lead- 
ing to stationary patterns will appear, (b) for q < and 
i] < (1Z ) only bifurcations with soft mode instabil- 
ities can occur, (c) for q < and 77 > (1Z |_) hard- 
mode instabilities with or without temporal oscillations 

may be found and (d) for q > and r) < (7£_| ) the 

homogeneous, stationary state remains linearly stable. 
In the fast reaction limit, where the time scale T _1 be- 
comes much shorter than other time scales (except the 
RC-relaxation time of transmembrane potential fluctua- 
tions), only soft-mode instabilities remain. By varying 
the ratio D = D c /D of diffusion constants or q — q c /q 
of effective electrophoretic charges within physiologically 
plausible bounds, it is possible to switch between soft- 
and hard-mode instability in and between hard- 
mode and oscillatory instability in 1Z |_. 

For soft-mode instabilities, we have derived an equa- 
tion of motion of slow modes near the transition, which 
is of the form of a dynamical Cahn-Hilliard equation. 
The transition will generically be discontinuous (back- 
ward bifurcation) and the non-linear evolution will be 
characterized by regions of nucleation and of spinodal de- 
composition. Quenches into the supercritical regime are 
predicted to show coarsening behavior with a coarsening 
length scale growing as i 1 / 3 . For hard-mode instabilities, 
an amplitude equation is obtained as in O and parame- 



APPENDIX A: GRADIENT EXPANSION NEAR 
SOFT-MODE INSTABILITIES 

We transform (^H) to the new variables 

u = (7+0 + 7-C)/r 

v=(g-0/r (Al) 

and get 

Td t u= (7 + + L>7_)V 2 7i + 7+7 _(f-L>)V 2 ?j + 
+ (7+ + r7_)V 2 $ + (7+ + r7_)V(wV$) 

+ 7+7-V(tA7$) (A2) 

as equation of motion for u. Note that all terms on the 
rhs of the equation for u are at least 0(V 2 ) as a direct 
consequence of channel number conservation. For v we 
get 

Td t v = -T 2 v + (7- + L»7 + )V 2 w + (f - D)V 2 u + 
+(1 - r)V 2 $ + (f - r)V(uV$) 

+ (7_ + r7 + )V(t;V$) (A3) 

Note that v decays with a rate O(V ) due to the —T 2 v 
term. Thus v is fast near a soft-mode transition and may 
be adiabatically eliminated. It does not enter the leading 
order of the gr adi ent expansion, which is obtained by 
replacing $ in ( |A^ ) from the leading (O(V )) order of 

= (1 - V 2 )$ + r)(u + 7_?j) + a(u + 7_u)$ . (A4) 

This leads to 

d t u = t7 x V 2 jr7 c w+ (I + u)cj){u) - J (j)(u)du)\ (A5) 

with t" 1 = (7+ + £>7_)/T and </3(u) = —iju/(l + au). 

The linear 0(V 4 ) term is obtained from the expansion 
of Lu+i^ k 2 ) = T -i{P( v ~7 h )/r lc + e s k 2 + 0(k*)}. For^ 2 
we obtain 

f = ~I + IF [ (1 " D ~ 11)2 + 2ll{l - ~ 7+) ~ A ^- Dq \ 

-^lh--J+)(l-D-T 1 ) + 2r 11 _Dq] 2 (A6) 

Adding a term — V u on the rhs of (|A§) we get @ 
and (§(]). 
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APPENDIX B: AMPLITUDE EQUATION NEAR 
HARD-MODE INSTABILITIES 

The system of nonlinear partial differential equations 
([|||) is written as £u = J\f(u, u) In the vicinity of rj c , k c 
u = (g, C, $) contains both fast scale variations (with 
Vit = 0(k c )) and slow scale variations (with Vu — > for 
V Vc)- The slowly varying part may be obtained from 
a perturbation expansion in e = y/\ (77 — r) c )/r} c \, which 
explicitly separates slow scale variations from fast scales 
by introducing appropriately scaled length and time vari- 
ables X = ex T = e 2 t. For one-dimensional patterns, the 
direction perpendicular to the pattern wave vector should 
be scaled as Y — e x l 2 y |18| ]. As functions now may de- 
pend upon fast and slow variables, we have to replace 
dt — ^ dt + £ 2 9t and d x — > d x + edx etc. In this way, 
the linear part £ becomes £ = £0 + e£\ + e 2 £ 2 - Insert- 
ing the expansion u = eui + e 2 u 2 + e 3 u 3 + 0(e 4 ) (with 
ui ( X, Y, T, x) — (x(X,Y,T)exp(ik c x) + cc) and sorting 
with respect to powers of e one gets a hierarchy of equa- 
tions, which allows to determine u„ from the with 
k < n. It starts with £o u i — 0, which implies that Ui is 
an eigenstate to £q with eigenvalue zero. As dr carries 
an explicit factor of e 2 , the equation, which contains the 
slow time derivative of Ui appears in third order in e and 
has the form £ 2 Ui + £iu 2 + £o u 3 = A/^Ui, u 2 ). Insert- 
ing U2 as obtained from the perturbation expansion in 
terms of ui, one gets the amplitude equation by taking 
the scalar product of the third order equation with the 
left eigenstate of £q to eigenvalue zero, u}. This solvabil- 
ity condition gives 

J u t 1 £ 2 uid 2 r + J u t 1 £ 1 u 2 d 2 r = J u\N z d 2 r (Bl) 



The Ihs of (Bl) is linear in \ and can be obtained directly 
from expanding uj + (rj.k 2 ) around r) c ,k c , io+(rj,k 2 ) = 
T^ h 1 {e 2 + (f l (k — k c ) 2 + 0[{k~k c ) A ]}. Inserting this expan- 
sion on the Ihs of (Bl) and neglecting the y dependence 
for simplicity gives 

r h d t - ^ - e h dl) X = J u|A/- 3 d 2 r (B2) 



For the time constant 77, we get 

T + (l + D)k 2 - Vc k 2 /(l + k 2 ) 



Th 



Dk*(k* + <y-+<y + /D) 



(B3) 



and the square £ 2 of the correlation length is 



u- 2 



(1 + k 2 ) (k 2 



+ 7- + ^) 



Evaluating the non-linear term on the rhs of (B2) gives 
(E2T) with the following expression for J 



J = 



(fc2+ 7 _ + 2)(l + fc 2 ) 2 in tn '- + D 
■{pi +<$> k ni + 2a) 



(l + k 2 )[2 Pl (?± + k 2 c +C 



-$ fec { ni { 1+ + k 2 ) - <n_(i + k 2 c ) Zl } 



The expressions ni, Z\ and pi stem from 0(e 2 ) correc- 
tions and have the form 



1 



7ll = 



18fc4 



{l + Akl){ 1+ /D + ik 2 c + Ck c q-/-) 



Zl 



-2a( 7+ + (?7 _+4fc 2 )|, (B5) 
m((3 + 2fc 2 )(l + 4fc 2 ) + 4fc 2 (a - mr, c ) 



Pi = 



/3(l + 4fc2)(l + fc2) 
fl + k 2 ) Em + a 



1 + 4fc2 



(B6) 
(B7) 
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4 l, 
FIG. 1. Different biological models. Left: simplified 

synapse. Right: Model of an axon. 
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k C) 

FIG. 2. Neutral curves in regions 1Z++ (a), and 1Z~- (b) 
for different values of q. For the upper, middle and lower 
curve q = 0.1, 0.2 and 0.6 in (a), q = -3.0, -2.0 and -1.0 
in (b) and q = —0.3, —0.5, —1.0. Other parameters are 
7+ = 0.01, 7_ = 0.5 and D = 0.5 
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FIG. 3. Boundaries between hard-mode and soft-mode in- 
stabilities in D — 7-plane in 1Z++ for different values of q. In 
regions below each curve k% 0. From lower to upper curves, 
values of q are 0.0, 0.1, 0.25 and 0.5 from lower to upper 
curves. 
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q = -2.0, 7 + = 0.01, 7 - = 0.5. P b) 



FIG. 5. The critical frequency £l c : a) As function of D for 
r = 4.0, 4.5, 5.0 and b) as a function of r for D = 2.75. The 
remaining parameters are 7+ — 7- = 1 and q = — 2. 



9 







-0.005 

flu) 









1 



















-1.5 -1 -0.5 0.5 1 

u 

FIG. 6. Local part of the effective potential f(u), which 
controls the Cahn-Hilliard dynamics near soft-mode instabil- 
ities at a = 0.6. The uppermost curve displays the situation 
at the transition (r) c /rj = 1.02525), whereas the middle curve 
corresponds to the spinodal point rj c /rj = 1. The lower curve 
is in regime at rj/rj c = 0.97. For better visibility, the f(u) 
at the transition point has been enlarged by a factor of 10. 
If a < 0.5 the global minimum of the potential is attained 
for positive u and and the graphs look like mirrored at the 
u — 0-axis. 




D c) D d) 



FIG. 7. J = 0-lines within the fcj? > 0-domain for g — 0.1, 
7 = 2.0 and various a: (a) 0.25, (b) 0.5, (c) 0.6 and (d) 0.75. 
The value of a crucially determines the number of domains 
with continuous transitions (shaded). If a > 0.5 the system 
clearly displays reentrant behavior. 
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